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INTRODUCTION 

Developments  of  modern  science  have  confronted  Che  scientist  and  the 
engineer  with  a  variety  of  problems  which  cannot  be  solved  formally.   Hence 
the  widenning  interest  in  numerical  analysis,  a  branch  of  mathematics  which 
leads  to  approximate  solutions  by  repeated  application  of  the  four  basic 
operations  of  algebra.  Numerical  analysis  has  been  applied  to  scientific 
and  technological  problems  from  the  very  beginning  of  applied  science,  but 
has  been  given  new  impetus  by  development  of  the  electronic  digital  computer. 
Many  classical  numerical  analysis  techniques  are  available  for  engineering 
applications  to  solve  differential  equations.  There  is  another  group  of 
procedures  unique  to  the  engineering  literature  which  generates  approximating 
recurrence  relations  from  approximate  Z-transforms.   No  matter  how  different 
all  these  techniques  are  in  terms  of  their  mathematical  approaches  and  the 
algebra  involved,  they  all  have  two  things  in  common;  the  calculations  are 
performed  with  discrete  values  and  on  a  step-by-step  basis.   Consequently 
the  time  interval  between  steps  is  an  Important  factor  which  affects  the 
accuracy  as  well  as  the  speed  of  the  computation. 

Up  to  a  point,  the  smaller  the  step  size  used,  the  longer  is  the  solution 
time  required  and  the  more  accurate  is  the  solution  obtained.  Since  signifi- 
cant digits  are  limited  in  computation,  the  increased  solution  accuracy 
produced  by  the  reduced  step  size  is  lost.  This  effect  is  termed,  "round- 
off error"  Therefore  there  exists  an  optional  step  size  to  be  used  in 
producing  the  numerical  solution  of  optimal  accuracy.   In  general  no  single 
technique  is  best  in  all  cases.   In  fact,  the  effectiveness  of  approximate 
methods  hinges  on  the  type  of  function  in  question  and  the  goodness  of  each 


method  Is  measured  by  a  number  of  factors,  namely,  the  program  set-up  time, 
the  solution  speed,  and  solution  accuracy  and  stability.   In  other  words, 
it  is  a  consideration  of  economy  and  accuracy  and  these  two  quantities  are 
contradictory  in  nature.  Therefore  a  brief  examination  of  strong  points  and 
weaknesses  of  various  types  of  digital  techniques  and  comparisons  among  them 
would  provide  a  guidepost  for  selection  of  an  optimal  technique. 

The  objective  of  this  report  is  to  review  previous  methods  for  digital 
.ution  of  differential  equations,  and  to  illustrate  their  applications  for 
approximating  the  solutions  of  ordinary  differential  equations. 


II.   Survey  of  tethods  for  Digital  Solution  of  Differential  Equations. 

Methods  for  approximating  solution  of  ordinary  differential  equations 
are  based  on  the  principle  of  discretization.   These  methods  have  the  common 
feature  that  no  attempt  is  made  to  approximate  the  exact  solution  y(t)  over 
a  continuous  range  of  the  independent  variable.  Approximating  values  are 
sought  only  on  a  set  of  discrete  points  t  ,  t.,  t2»  ,,,tn«   Generally  speak- 
ing, a  discrete  variable  method  for  solving  a  differential  equation  consists 
of  an  algorithm  which  furnishes  a  number  y  corresponding  to  each  lattice 
point  tn.   Which  is  to  be  regarded  as  an  approximation  to  the  value  y(t  )  of 
the  exact  solution  at  the  point  t  . 

Discrete  variable  methods  fall  into  two  classes:  classical  numerical 
analysis  techniques;  and  numerical  transform  techniques.   Classical  numerical 
analysis  techniques  yield  one-step  methods  and  multiple  step  methods.   In 
a  one-step  method  the  value  of  y_  can  be  found  if  only  one  initial  value  is 

It 

known.   In  a  multiple  step  method  the  calculation  of  y  -  requires  explicit 
knowledge  of  more  than  one  starting  value.   Oftentimes,  a  one-step  method  is 
used  to  start  a  multiple  step  method. 

The  numerical  transform  techniques  were  first  introduced  by  Tustin  and 
were  further  expanded  by  Kadwed,  Boxer-Thaler,  and  Halijak.   These  techniques 
generate  approximating  recurrence  relations  from  approximate  z-transform  of 
l/sn  and  f/sn. 

III.  Classical  Numerical  Analysis  Techniques. 
1).   One -Step  Methods 


Lat  a  typical  first-order  differential  equation  be  given  by 


dy 

-—  -  f(t,y) 

dt 


y  -  yo     at  t  -  tQ  (3.1) 


and  let  it  satisfy  Lipschitz  conditions  in  soma  closed  region  D.  There 
exists  a  single-valued  function  y(t)  continuous  in  D  such  that  it  satisfies 
Eq.  (3.1),  then  it  can  be  solved  approximately  by  these  methods. 

One-step  methods  are  usually  divided  into  two  classes;  The  first  class 
includes  Taylor  series  method  and  Picard*s  method.   In  these  methods,  y  in 
Eq.  (3.1)  is  approximated  by  a  truncation  series,  the  individual  terms  of 
which  are  functions  of  the  independent  variable  t.  The  second  class  is  re- 
presented by  the  methods  of  Euler  and  Runge-Kutta  methods. 

A.  Taylor  Series  Kathod. 

Consider  the  differential  equation  (3.1)  with  the  initial  condition  y-y 
at  t»t  .  Lat  the  required  solution  be 

y  -  y(t)  (3.2) 

If  t»t  is  not  a  singular  point  of  the  function,  y(t)  can  be  expanded  in 
Taylor  series  about  this  point.  Thus  with 

(m) 


o 


yo 


d  y 


dtm 


(3.3) 


C )   ^       2  (2)   ^        (3) 
y  -  yQ  ♦  (t-tQ)y;   +  ~<t-t0)  yQ   +  -^oK        +     "  (3'4> 


a  power  saries  in  t  that  converges  over  some  range  t  <t<t  . 

o      o 

The  value  of  y    is  evaluated  by  making  use  of  the  differential 


equation  (3.1)  which  can  be  written  as 


y(1)  -  f(t,y)  -  gl(t,y),  (3.5) 


■c-n-ci-jing  Eq.  (3.5)  yields 


(2)  2&i     3gL   /,%  q\ 

y   -  — ■  ♦  —  yu;  -  s2(t,y,yu;),  (3.6) 

at  dv 


and  by  repeated  differentiation 


where 


y<ra>  -  gm(t,y,y<1>,y<2)>...y<m-1>)         m  -  1,2-.      (3.7) 


'm-1   .  m:2   agn,-l    (1+1) 


^t       i-0   ^y' 


Setting   t-tQ  and  y-y  ,  Eq.  (3.7)  becomes 


"■VVV^.-','-11'.         -».«•••  o.9) 


The  truncation  error  after  the  mth  terra  is  given  by 


R 


m 


fm<?) 


m  t 


m 
(t  -  to)  ,     where     tQ<  ^  < 


(3.10) 


Replacing  f (£)  with  an  upper  bound  M    to  its  value  in  the  interval 

(t,  t  ).   The  truncation  error  is  then  bounded  by 
o 


R  < 

ra  — 


M 


(ra) 


ra 


<t  -  to> 


in 


(3.11) 


In  particular,  for  a  convergent  Taylor  series  with  alternating  sign,  the 
truncation  error  after  rath  term  cannot  exceed  the  (m«s-l)th  term  Eq.(3.11) 
becomes 


R  I  < 
ra 


f     (t  ) 


(t  -  t0) 


m 


(3.12) 


B.   1.   Euler's  Method. 


This  method  is  of  very  little  practical  importance,  but  it  illustrates 
in  simple  form  the  basic  idea  of  those  numerical  methods  which  seek  to 
determine  the  change  of  A  y  in  y  corresponding  to  a  small  increment  of  the 
argument. 

Consider  Eq.(3.1),  the  left-hand  side  of  the  first  part  is,  by 
definition 


dy           A.y 
-   lim  (3.13) 


dC     ^t-»0  At 


therefore,  for  small  value  of   t 


dy 
Ay   -  -<it 

dt 


Thus,  the   increase   in  y(t)  when  t   increases  to  t+At   is  approximately 


ym.1  -  y„  -  f<t  »  ym)T  (3.14) 

...••..         in  mm 


where 


T  -  Zi  t 


However,  the  exact  expression  for  y  at  t-t     ,,  using  Taylor  series  formula 


z       -z     +  Tf  (t   ,  a   )   +  -  T2y<2)(J)     .  (3.15) 

m+1     m  o  m+l 


\.:.~scj 


'■      <   $     i  <  t 
m        >m+l        m*i 


m-J-1 


Subtracting  Eq.(3.14)  from  Eq.  (3.15)  results  in 


z  ,-s  -  y  ^  -  <z    -  y  )  +  T<f  <t  ,  z  )  -  f  <c  ,  y  )) 

m+i         m+1  mm  mm  mm 


+  -  T     y     \%      V  (3.16) 

2  m+1 


The  left-hand  side  of  Eq.(3.16)  Is  by  definition  the  truncation  error  of  y 
at  t»t 


'm+1 ' 


Let  6  -  z  -  y  (3.17) 

m+1  m+1  m+1 


f(t      ,  z    )   -  f(t    ,  y  )        ^f(t    ,  z    ) 

mm                mmw             mm  /•>    io\ 

and  k     ■  »  Vj.IoJ 

m  z     -  y  °  y 

mm 


Eq. (3.16)   becomes 


1       o        /o\ 

6     _  -e    +  Tk  ^    +-t    y      (H  ^).  (3.19) 

m+1  m  m  m       2         *  m+1 


For  m  -  0 


t0  -  0  (3.20) 


m  »  1 


^x  -  -  R2  y(2)(^),  V  ^<  tx  (3.21) 


a  =  2 


1     o      (2)  ^     9      o\ 

e2  -   (t   +  TK1)-  T     yV        (^)    +  -  TZ  yU;    <%)  (3.22) 


where 


> 


tl<f2<t2      J 


proceeding   in  this  manner,  the  truncation  error  at  nth  step  is 


where 


£      -  T2 
n 


,1) 


n-1 


(2) 


;.-: 


(£>     tt     (1   +  TK  )   +  -  yv""    (f2)     tt     (1   +  TK  ) 

jol  J  9  ^«9  J 


J-2 


+   ...     +  -   y(2>     (fn) 
2 


(3.23) 


or 


n-1 
£_   -  T-       Z       A,   a, 

n  :-G        X      i 


(3.24) 


n-1 
A.    «       TT      (1   *  TK.) 

x     i  J 


A,-l 


i  *  n 


i   -  n 


-.'.a 


a     -  -  y 
2 


(2)  <n> 


(3.25) 
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Thus  errors  -  yn(|:i)  introduced  at  each  step  because  of  the  inaccuracy  of 
formula  (3.14)  are  to  be  multiplied  by  amplification  factors  A.  before  being 
summed. 

C.  Runge-Kutta  Mathod. 

This  method  is  an  algorithm  designed  to  approximate  the  Taylor's  series 
solution. 

Consider  the  Taylor's  Series 


T2 

r    .,  -  y  +  T  y  +  -  y  +  •••,  (3.26) 

n+l    n     n   «   n 


where 


y  .,  -  y((n+l)T)     and     y  -  y(nT). 
n+l  n 


If  the  Runge-Kutta  formula  is  derived  by  reta inning  terms  in  the  Taylor's 
series  expression  up  to  mth  power  of  T,  this  formula  is  called  the  Runge- 
Kutta  method  of  mth  order  accuracy  or  the  Runge-Kutta  mth  order  method. 

C.   1.  Second  Order  Runge-Kutta  Formula, 

Consider  increments  in  y  defined  by  the  equations 


^y  -  R^'y  +  R2A"y  (3.27) 


where 


^'y  -  f <tn,  yn)  T, 


11 


A"y  -  f  (tn  +  aT,  yn  +  bT)  T, 


and      R.  ,  R~>  a  and  b  are  constants. 

Expanding  Eq#(3.27)  with  respect  to  Taylor's  series  of  two  variables, 
results  in 


Ay  -  T(R1  +  R^  fn  +  (a(fc)n  +  b(f  )R)  TZ   R2  +  •• 


(3.28) 


where 


and 


fn-f(Cn'yn>  "  *n  • 


2t 
(0„  -  (— ) 


tatn*  y-v 


3f 

(f  )   -  (— ) 

y  n      ' 

3y 


tmtn>  y-yn 


Eq.  (3.28)  will  be  equal  to  Eq.(3.26)  up  to  second  order  of  T  if  its 
parameters  take  the  values 


Rx  +  R2  -  1 


(3.29) 


R2a  -1/2 


R2b  -  1/2 


There  are  four  constants  to  be  determined  and  only  three  equations.   Choose 
R2  as  parameter,  the  constants  a,  b,  and  R,  can  be  determined  in  terms  of 
R2.  Thus 


12 


a  -  b  -  (3.30) 

2R2 


Rx  -  1  -  R2 


Consequently,  an  infinite  number  of  forms  of  Eq. (3.28)  can  be  established. 
Formulas  of  higher  order  can  be  obtained  in  the  similar  manner  although  no 
formula  beyond  fifth  order  has  ever  been  developed. 

C.   II.  Third  Order  Runge-Kutta  Formula. 

Again  take  the  differential  equation 


dy 

»  f(t,  y).  (3.31) 

dt 


and  expand  it  with  respect  to  (t  ,  y  )  by  Taylor  series  for  two  variables 

n   n 

to  obtain 


ay  -  fnT .+  -  Cft  ♦  fyf)n  t2  *  -  (ftt  ♦  2fcy  f  +  fyy  f: 


f3 

t    *y  "'  *y 


')n  : 


Define  the  increment  of  y  by 


Ay  -  R^'y  +  R2^"y  +  R3^*"y  (3.33) 
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where 


^'y  -  f <tn»  yn)  T 


Any  -  f(tn  +  iaT,  yn  +  mzQ»y)   T 


A*"y  -  fCt,,   +  XT,  y^  +  p^"y  +  (X  -  p)  *  »y)   T 
n  o 


The  symbols  m,  A  >  £nd   p  are  constants.      Expanding  Eq. (3.33)  with  respect   to 
Taylor's  series  at    (t   ,  y_)»  results   in 


*'y  -  Tfn 


A-y  -  T{fn  +  nT(ft)n  +  B4'y(fy)n  +  —  [(mT)2(ftt)n   +  2m2T^'y(fty)n 


+  M«y)2(fyy)n  ]  + 


A'"y  -  T^  fn  +^T(ffc)n  +  [/^"y  +  (X  -/*)  A<y  j  <f  )n 


1 
2 


-{ aT)2(ftt)n  +  2xt[/>  ^»y  ♦  a  - n  A'y  j  <Vn 

+  [f*"y  +  a  -/°)^,y  ]2<fyy)n)  +  •••  )         <3-34> 


where 


32f  S>f 

<?y     u  cn,  y  yn#  y  ay  t-tn,  y«yn 


5>f 
(ft)      -    <— )  .  f      -  f(t    ,   y    ) 

at  t«tn,   y»yn  n  n       n 
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Eliminating  ^ *y  and  zJ'y   from  the  right  hand  side  of  the  second  and  third 
equations  of  Eq. (3.34)  and  substiting  into  Eq.(3.33),  the  result  coincides 
with  Taylor's  series  expansion  up  to  the  third  power  of  T  provided  the 
following  relations  hold  among  the  constants 

Rx  +  Rg  +  1*3  -  I  (3.35) 

1 

R~m  +  R3X  a  - 
2 

2       9    1 

R  mz  +  R  \*   *>   - 

*  3 


1 

>m  si 

3> 


Rofm 


Since  there  are  six  constants  to  be  determined,  namely,  R  ,  R  ,  R  ,  m,  \9 

12   3 

and  f  9   and  only  four  equations,  the  constants  m  and  X  can  be  chosen  as 
parameters.   The  constants  R-,  R2,  R,  and  p  as  calculated  from  Eq.(3.35)  are 


6m\  -  3(m  +  \)   +  2 
Rt   -  (3.36) 

6m\ 


2  -  3X 


R 


2 

6m(m  -  \) 


2  -  3m 


R3 


P 


6\(\   -  m) 

\(\   -  m) 
m(2  -  m) 
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From  Eq.  (3.35),  it  is  clear  that  R~  cannot  be  zero  and  R.  and  R2  cannot 
both  zero.   If  R2«0,  then  by  Eq.  (3.35) 


1  2 


\   -  -     IL  -  -     R.  -  -     and  fm   ■  -  (3.37) 

3      J   4      X   4  3 


Itoreover,  assume    f  to  be   1/3,  then  m  is  equal  to  2/3.     Therefore 


1 

Ay  »  -  (^*y  +  3^«"y)  (3.38) 

4 


^•y  -  f(t   ,  y  )  T 

n       n 


^"y  °  f  (t     +  -  T,  y     +  -   A*y)  T 
3  n       3 


2  2 

At..y  -f  [t     +-I,  y     +.(A'y+  ^»y)    ] 


On  the  other  hand,  if  \«0,  then 

2  3        1  y>-  1 

m  =  -     R  -  -     R«  —     and     R.  -  (3.39) 

3  Z       4      3   4f  V 

If  it  is  assumed  that  f-lt   then  another  form  can  be  constructed  from 
Eq.  (3.39)  and  Eq.(3.33);  it  is 


1 
Ay   -  -  O  A"y   +  ^»"y)  (3.40) 

4 


where 


^•y  -  f (tn,  yn)  T 
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A,,y«f(t  +  -  T,  y  +  -  A  *y)  T 


n   3     n   3 


A»»y  -  f  (t  ,  y^  •>  a  "y  -  a.  «y)  T 


C.   IV,  Fourth  Order  Runge-Kutta  Formula. 


The  Runge-Kutta  fourth  order  formula  has  the  form 


A  y  -  Rj  2i9y  +  R2^«8y  +  R^^'y  +  R4^m'y  (3.41) 


where 


<a'y  -  f  <tn,  yn)  T 


A"y  »  f  (t  *  mT,  y  +  m  -d*y)  T 


^«»y  «  f(tn  +  XT,  yn  +  /^"yn  +  (X  -/?)^cy)  T 

^imy  „  f(t   ^Ts  y   +^^«iy  +^^»y  +  (^/-  (^-^)Zi8y)  T. 
Expand  second,  third,  fourth  and  fifth  equations  of  Eq. (3.41)  to  obtain 


A'y   -  f  (t  ,  y  )  T  (3.42) 

n   n 
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A»y  -  T(f(tn,  yn)   +  «nT(ft)n  +  m^«y<fy>n 


+  —  ((ml)2    (ft|;)n  +  2(mT)<ra  A«y)(fty)n  +  (xn  ^»y)2(fyy)n) 


+  —  f(mT)3    (f_^)      +  3(mT)2(m^,y)(ft.4.   ) 
«l    l  ttt  n  J       tty  n 


+  3(EiT)(m^«y)2(f       )      +  (n^«y)3(f       )     1  +  •••  ) 

cyy  n  yyy  n  j  ) 


A*"y  -  T  (f(tn,  yn>  +  XT(ft)n  *  (/^"y  +  a  -/>)2i,y)<fy)n 


+  _  [aT)2(ftc)n  +  2aT)<f^"y  +  a  -/°)  ^y><fty>n 


*  <fz>"y  +  a  -P)^»y)2(f     )    ] 

'  '  yy  nJ 


+       j"  aT)3(fttt)n  +3aT)2(/^«y  +  a  -/O^yXf,.^ 


+  3aT)(/^"y  +  a  -/°)^'y>2<ftyy>n  +  o^"y  +  &  -z^'y)3 


+  (/^"y  *  (X  -/O^'y)   (fyyy) 


J—} 


^<4)y  -  T  (fn  +/vT(f,c)n  +  (<?~A(3>y  +  T^i2\  +  {// -  ^-^O^'yXf  >n 


+  —  [^T)2(ftt)n  +  2^T)(^3)y  +  r^\<2)y 


w). 


+  (/^-  ^-^)^'y)(fty)n  +(<//-  <^  -^)Zfy  +  <?^v  'y 


^^(2)y)2(fyy)nJ  +  ... 
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Eliminating  ^9y9    ^"y  and    ^"»y  from  the  right-hand  side  of  Eq.  (3.42)  by 
successive   substitution,   putting  these  results   into  Eq. (3.41)  and  comparing 
with  Taylor's  series  yield 

Rx   +  R2   *  R3   *  R4  «  1  (3.43) 

1 
R  m  +  RA   *  RJA  *  - 

1 

R_m2   *  RA2   •>  RJ?  -  - 

2  3  *  3 

1 

R„m  f  ♦  R,  (mcA*  X  )   »  - 

34  6 

o  o  *5 

R  nT    *  R^X      *  R  /f*   -  - 
1  J  4  4 


RrnXf-i-  r  ^(bkt+  X£)  »  - 
J  4  8 


1 

R^m3f  +  R   (nj<r-+  X2*?)   -  — 

->  *     *  12 


1 
4  24 


There   is  a  one -parameter  family  of  solutions  derived  by  Kutta  as  follows: 


1  2  -  t  t  1 

R    -  -  R     "  R,  -  -  RA  -  -  (3.44) 

1       6  2  3  J       3  4       6 


m 


1/2,     X  "  1/2,     /«  l/2t,     ^»  1,     <?-  1  -  t,      T>  t 
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If  t-1,  a  Runge-Kutta  fourth  order  formula  can  be  constructed  as  bellows: 

Ay   -  -  <^'y  +  2^»y  *  2  A(3)y  +  ^(4)y)  (3.45) 

6 

where 

n 


^«y  .  T  f(t  ,  y  ), 

n   n 


A*  'y  -  T  f  (t  +  -  T,  y  +  -  A'y) 
2       2 


/3v  l  1 

^v  y  =  T  f(t„  +  -  T9  y  +  -  ^f:y) 
2     n   2 


(4) 

^v  :  -  T  f  <t  +  T,  y  *  a  »y) 
n      n 


4).  Fifth  Order, Runge-Kutta  Formula 

The  5th  order  Runge-Kutta  formulas  were  derived  recently  by  H.  A.  Luther 
and  K.  P.  Konen.  The  derivations  are  as  follows: 


I-t       dy/dt  »  f (t«  y)   together  with  y(t  )  »  y  (3.46) 

o     o 


The  fifth  order  Runge-Kutta  formulas  are  phrased  as 


t  (i) 

A  y  «■  +  2  R.  A      y 

1 
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A'y   -  T  f (tn>  yn) 


(3.47) 


,(s). 


s-1 


0), 


^v°'y  -  T  f  (tn  +  agT9  y  *  Z    b   ^vjyy). 


where  2<s<6.  and  the  coefficients  R.  are  the  constants  to  be  determined. 
The  usual  procedure  yields  44  equations  involving  the  various  parameters. 
Assume  that 


s-1 

a  "  2  b  ,, 

■   j-1  s^ 


2«s  ^  6 


(3.48) 


and  eliminate  easily  identifiable  combinations  to  obtain  the  following  16 
relations 


Rx   *  R2  *  R3  +  R4  +  R5  +  R6  -  19 


a2R2   *  a3R3   *  a4R4  *  a5R5   +  a6R6 


a2R2   +  a3R3   +  34R4  +  a5R5   +  a6R6 


3  3  3  3  3 

a  R     +  a  R     +  a  R     +  a  R     +  a  R 
2  2         3  3         4  4         5  5  6  6 


2 
1 
3 

1 
4 


4  4  4  4  4  1 

aoRo   +  a<A    +  a  R     +  a  R     +  a  R     -  - 
22         33         44         55  66       5 


(3.49a) 
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CIR3   +  C2R4  +  C3R5   +C4R6  "7 

o 


a.c.R.    +  a.c.R.    +  a_c_R_    +  a.c.R.   »  1/8, 

3   i3  424         5  J  5  64o 


♦  A   *  V2R4  +  a5C3R5   +  a6c4R6  "  1/10« 


d.R^    +  d_R.    +  d_R_    +  d.R,   -  1/12, 
-^  24  35  46 


S3dlR3   *  V2R4  *  a5d3R5   +  VA  "  1/15' 


CA  +  SR/  +  C<A  +  c/ra  "  1/20» 

13  24  j>   5  46 


(3.49a) 


ab„   R     *  (a  b     +  a  b  0)R,    +   (a  b  rt    +  a  b  „    +  a  b   ,)   R- 
2  32  3  2  42       3  43     4  2  52         3  53         4  54       5 

3  3  3  3 

+   (a  b  ^   +  a  b  „    +  a  b   ,    +  a_b     )   R  -  1/20 
2  62         3  63         4  64         5   65 


°lb32R3   +  felb5S   +  C2b54)   R5   *  <Clb63   *  C2b64  +  °3b65>  R6  "  1/24' 


d.b.^R,    +(db        +db)R+(db        +db        +db)R-  1/60, 
j        1  43  4        v  l  53         2  54       5  1  63         2  64         3   65        6  * 


<a3   *  a4)clb43R4  +  [(a3   +  a5)clb53   +  (a4  +  WsJ 


54 


+     [(a3   +a6)cib63   +  (a4  +  a6)c2b64  *  (a5   *  a6)c3b65]R6  "  7/12° 

Clb43b54R5   +   [  Clb43b64  +  (°lb53   +  C2b54)b65 J  R6  "  1/12°>  J 

(3.49a) 
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where 


1+1  i+1       2 

c*    "       2       a<b<^      i    »  d.    -       S       a  b,    2      .   .  (3.49b) 

1         j-2       J        2*  J  j-2       J  '  J 


simplification  requires  that 


R2  -  0  (3.50) 


2 

ai 
and  ~»c,   „   ,  3$i^6  (3.51) 


'i-2    f 


then  Eq. (3.49)  can  be  simplified  considerably.  After  eliminating  duplicates 
and  combining  some  equations  solve ,  in  addition  to  (3.51), 


R,  +  R0  +  R,  +  Rc  -2-  R,  =  1,  (3.52a) 

13    4    5     6 


and 


a3R3   *  a4R4  *  a5R5   +  a6R6  "  1/2» 


2  2  2  2 

a0R0    +  a.R,    +  a^R,.    +  a,R„   ■  1/3, 
33  44  55  oo 


3  3  3  3 

a^   +  a4R4  +  a5R5   +  a&R6  »  1/4,  (3.57b) 


a3R3   *  a4R4  *  a5R5   *  4R6  "  i/5> 


end 
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and 


a3b43R4  +  (a3b53   +  a4b54)   R5   +   (a3b63   +  a^  +  a^)   R6  -  1/6 


a3b43R4  *  (a3b53   *  a4b54>  R5    *  (a3b63   +  a4b64  +  a5b65)  R6  "  l/12 


a3b43R4  *  (a3b53   *  a4b54}   R5   *  (a3b63   +  a4b64  *  a5b65>   R6  "  1/2°' 


V3b43R4  +  a5  (a3b53a4b54)   R5   +  a6(a3b63   *  a4b64  +  a5b65}   R6 

-  1/8, 

a4a3b43R4  *  a5(a4b53   *  a4b54)  R5   *  a6(a3b63   +  a4°64  *  a5b65}  R6 

-  1/15, 

<3.52c) 


a3b43b54R5  +  [a3b43b64  *  (a3b53  +  a4b54}  b65  J  R6  °  l/24> 

2  r  2  2      2 

a3b43b54R5  *  L a3b43b64  +  (a3b53  *  a4b54>  b65  J  R6  °  l/6°* 

(3.52d) 

The  situation  is  now  as  follows.  Equations  (3.52b),  (3.52c),  and 
(3.52d)  are  to  be  solved  independently.   Then  (3.52a)  yields  Ri#  Equation 
(3.51)  are  used  to  find  b32,  b,2,  b<-2»  ^52*  Then  equations  (3.48) 
determine  b_,,  b_,,  b. , ,  b__  and  b._.   This,  with  Ro«»0,  completes  the  solution. 

Zl    jl    <+l    jl       01  1 

The  family  of  solutions  due  to  Kutta  may  now  be  found  by  taking  b^c^O. 
From  (3.52d),  a  =2/5.   It  then  develops  that  the  third  equation  in  (3.52c) 
has  the  left  member  equal  to  -&3a4  times  the  left  member  of  the  first  of 
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this  group,  plus  a^  +a,  times  the  left  member  of  the  second  of  this  group. 
This  forces  a,  to  be  1.   Equations  (3,52c)  may  now  be  solved  for  h^/Rg  and 
b54R5*  Wnen   the  results  are  substituted  in  a^b^  ^54%  +  b64R6^  "  1/24 
(the  consequence  of  (3.52d)  and  h65°0)  it  is  found  that  b^  "IS/A.   In 
summary,  there  results  the  following  description  of  a  family  of  solutions: 

R2  "  °»    b65  "  °»    b43  "  15/4»    a3  "  2/5»    a4  "  ls»   (3*53> 

with 

R3  -  125  [lOa5a6  -  5(as   *  afi)   +3]  /  [  72(2  -  5a5>(2  -  5afi)  ]  , 

R4  »    [  8a5a6  -  7<a5   +  afi)   +  6  ]  /  [  36(1  -  a5)(l  -  afi)  ]  , 

R5  "    [  l   "  a6  ]  /  [12a5(1  "  a5)<5a5   -  2><a5   "  a6>]  »  (3.53) 

R6  "    [  l  '  a5]  '  [I2a6(1  -  a6)(5a6  "  2)(a6  "  a5}]  • 

Rl  "  l   *  R3  "  R4  "  R5   -  R6» 


cs.d 


b53   -  5  [7   -  10a6  -   108R4(1   -  a6)  ]  /  [  144R5(a5   -  a&)  ]  , 

b54°    I1   -*6M36R5<a5   -a6>]  • 

b63  "  5  [7   "  10a5   "  1°8R4<1  "  a5}J  fi  144R6(a6  "  V]  »        (3.53c) 

b64  "    f1   -a5M36R6<a6  -a5>]  » 


and 
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b32  "  2/  I  25a2  |  , 


b42  "  1/a2' 


b52  "  [  5a5  "  4t53  "  10b54]  '  ^  ^  ^  ^^ 

b6,  *  t5^-4^-10^]^^]. 


and 


b   -  a  -   2   b   ,         2^j<  6.  (3.53e) 


It  becomes  very  simple  to  construct  a  fifth  order  Runge-ICutta  formula 
based  on  the  coefficients  derived  froa  equation   (3.53).      The  result   is 


y  =  y     +(    4ABy  +   (16   +76)  A(5)y   +   (16   -  /6)  ^(6)y  )  /36, 

r.vl  n       I  J 

A'y  -  Tf(t   ,  y  ), 
n       n 

As,y  »  Tf  (tn  +  4T/11,  yn  +  4a  'y/ll),  (3.54) 


4'"y  -  Tf(tn  +  2T/5,  yn  +  |   9Vy  +  llA"y  }  /50), 

A(   ;y  -  Tf(t      +  T,  y     +   {   -llAs'y   +  15A,s,y   }/4). 
n  n 


A(5)y  „  Tf(tn  +  (6  -   /6)T/109  yn  +    {  (81    +  9s/6)  A*y 

+  (255   -  55v/6)z}'"y   +   (24  -   14/6)  A(4>y]/600), 

A(6'y  .  Tf(tn   +  (6   +  /6)T/10,  yn  +     {  (81   -  9^6)  A  'y 

+  (255   +  55^6)  a  «"y  +  (24   +  14/6)  z^4)y  "J  /6000) , 
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Let  b   =  0  instead  of  letting  b  _  «  0,  then   Eqs.  (3.5 2d) 


yields 


(a3b53  *  a4b54>  b65R6  "  1/24> 


<a3b53  *  a4b54)  b65R6  "  l/60> 


Using  these  equations  in  conjunction  with  the  first,  second,  fourth,  and 
fifth  equations  of  (3.52c),  the  terms  in  b63,  b64,  and  b65  catl  ba  e*ira*nated 
and  there  result  two  equations: 


R  (a6  -  a5>  -  (4a6  -  3)  ^^R^ 


R  (a  -  a  )  «»  (5a  -  4)  b  R  , 
v  6    5       6      65  6 


These   lead  to  a  >=1.     Another  family  of  solution  occurs  with 
6 


&0   -  0>  b7.    =0,  a     -  1,  (3.55a) 

i.  43  6 

and 

R     -   |    3  -  5(a4  +  a5)  +  iOa^  ']/   |60a3(a4  -  a3)(a3  -  a5)(a3  -  1)], 

R4  -  [  3  -  5(a3   +  a5)  +  10a3aJ    /  [  60a4(a3   -  a4)  (a4  -  a5)(a4  -  1)]  , 

R5   -  [  3  -  5(a3   +  a4)  +  lOa^  ]  /  [60a5(a3   -  a5)(a5  -  a4)(a5  -  1)  ]  , 

R     -  f   12   -  15  (a_    +  a.    +  a,.)    +  20(aoa.    +  a.a_    +  a.a.) 
6       L  3         4         5  34         45  53 

-  30a3a4a5  ]  /  [  60(1   -  a3)(l   -  a4) (1   -  a5>  ]  , 
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R,    -  1   -  R,    -  R,    -  R-    - 
1  3  4  5 


(3.55b) 


arid 


'53 


J54 
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[  2  -  5a4  ]  /  [  120R5a3(l  -  a,.)^   -  «4>  ]  . 
[  2  -  5a3  ]  /  [  120R5a4(l  -  a5>(a4  -  a^]  , 


6   -  2a,    -   10a,    -  14a_    +  5a_ 
L  3  4  5  3 


a,    +  25a, a_   +  10a 
4  5 


'64 


65 


-  20a5a4  J/[l20R6a3(l   -  ^^^3   -  a4)(a3   -  a^],      (3.55c) 

r  2 

6  -  2a.    -   10a„    -   14a-    +  5a_a,    *  25a0a_    +  10a. 

L  4  j  5  34  35  5 

~  20a5a3  ]'  [120R6a4(1  "  a5)(a4  "  a3)(a4  "  a5}]  ' 

[  3  -  5a3  -  5a4  +  lOa^  ]  /  [  60R6a5 (a5   -  a3)(a5   -  a4>  ]  , 


c\c 


'32 


42 


'52 


3  /  [  2a2  1  » 

4  '[2a2]' 


r  2 

L  a5   -  2a3b 


3  "35 


"  2%h5^fl2al\  • 


b62  "  t  l  "  2a3b63  *  2a4b64  "  2a5b65  J  '  i  2*2  1  • 


(3.55d) 


add 


j-l 

b   -  a  -   2   b  , 
Jl    J    i-2   jl 


2^j  ^6 


(3.55e) 


28 


Thus  i  t  can  be  constructed  as  follows: 


yn+l  "  Yn  5^8"y  +  5^(5)y  +  ^(6)y 


k,  ■  Tf  (t  ,  y  ), 
1      n'  Jn   * 

k2  «  Tf(tn  +  T/2,  yn  +  1^/2), 

k3  -  Tf(tn  +  (5  -  5)T/10,  yn  +  [  2k1  +  (3  -  5)k2  J/10), 

k4  -  Tf (tn  +  T/2,yn  +  kx  +  k2/4), 

k5  -  Tf  (tn   +  (5  +  /F)T/10,  yn  +  f  (1  -  1/5)k1  -  4k2(5  -J-  375)k3 

+  81«4  ]  /20, 
k6  »  Tf  <tn  ->  T,  yn  +  (  5  -  l)k1  *  (2/J  -  2)k2  -J-  (5  -  /5)k3 

-  8k4  *  (10  -  2/5)k5/4. 

This  fifth  order  formula  is  not  in  Kutta's  family.   It  belongs  to 
Lobatto  quadrature  formulas  which  have  errors  of  order  T  rather  than  T  . 

The  accumulated  truncation  error  of  Runge-Kutta's  method  can  be 
calculated  as  follows: 


Ijjc 


Zn+1  "Zn  *^z<tn,  \   J  T)  (3.57) 


be  the  exact  solution  of  Eq,  (3.46), 
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and  y         -  y     +  T*y(t   ,  y     ;   T)  (3.58) 

n*l         n  n       n 


be   the  calculated  value   of   the   solution 

Subowi- citing  licuj-c.cn    (3.57)   :rv;.  ^u^;.c:.    w„S3)   yields 


e 
n 


l+1  °*  en  *  T  i^(tn>  ^n   >   T)   '  <*z(V  an   '>   T)i  (3'59) 


by  application  of  the  triangle  inequality,  Equation  (3,59)  becomes 


en+li^    I   Snl     +TIUy<V  yn   ;   T)   "  * 2  <  V  2n   J   T)  H 


+  T   |Uz(tn,  zn   ;   T)   -  Aa(tn,  zn   ;   T)  ||  (3.60) 


The   Lipschitz  condition  yields 


^y(t   ,  y     ;  T)   -  4z(t   ,  z      ;  T)|l^   Lily     -  z 

n'     n   '  n'     n   '        U  II    n         n 


-L||e     ||  (3.61) 


and  |j^z<V  zn  ;   T)   -  ^s(tn,  zn   ;   T)|U  N(TP)  (3.62) 


where 


(p) 
N  -  Max   ||  f        (y) 

(P  +  1) ! 


and  p  is  the  order  of   the   Runge-Kutta   formula 
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Hence  from  Equation  (3.60)  the  required  truncation  error  is 


•»♦!  II  ^U  *LT)  'I  en  H  +  T**"  I-  0.63) 


The  value  of  L  can  be  determined  from  a  Runge-Kutta  formula  and  different 
formulas  correspond  to  different  values  of  L.  All  L's  have  upper  bounds. 
The  Runge-Kutta  method  seems  to  be  tedious  because  values  of  f(t,y) 
have  to  be  calculated  a  number  of  times  per  time  increment.   However 9  the 
formulas  are  systematic  and  hence  can  be  easily  programmed  on  an  automatic 
machine.  No  special  starting  procedure  is  required  and  calculations  can 
often  be  checked  by  repetition  using  a  different  step  size.   Furthermore, 
such  a  method  is  particularly  useful  if  certain  coefficients  in  the  differen- 
tial equation  are  empirical  formulas  for  which  analytical  expressions  are 
not  known.  The  step  sis©  can  be  altered  as  desired.   It  should  be  understood 
that  the  derivatives  as  evaluated  by  Runge-Kutta  process  are  not  the  actual 
derivatives  at  the  various  points  within  the  step  as  commonly  assumed. 

2).   Multiple  step  Methods. 

The  one-step  methods  are  necessary  for  obtainning  initial  values  in 
the  solution  of  a  differential  equation.   However9  they  involve  too  much 
labor  to  be  used  for  obtaining  a  numerical  solution  over  an  extended  range. 
This  can  be  offset  by  using  multiple  step  methods. 

Multiple  step  methods  are  always  expressed  in  the  difference-differential 
equation  form 
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°k  W  *  afc-I  ^--:   -•••    *  %  yn  -  T  [h  fn+k  +  *k-l  fn,:cl 


+  •  •  •    +  3     f 
o     n 


n  -  0,   1,  2,   •»•  (3.64) 

where  k  is  a  fixed  integer,  f  «■  f(t  ,  y  )(m  »  0,1,3,»««),  and  where  expand 
g  (//»  0,1,2, •••)  are  real  constants  which  do  not  depend  on  n.  Any  a.  is 
always  assumed  unequal  to  zero.  Equation  (3,64)  defines  the  general  linear 
k-step  method.   If  B^O,  the  formula  (3,64)  is  called  "closed";  otherwise  it 
is  called  open. 

Unlike  one-step  methods,  multiple  step  methods  are  not  self -starting;  if 
some  values  y-^ijy.,.  o>,,#y  are  r*ot  known*  these  methods  break  down. 
Such  is  the  case  at  the  beginning  of  the  computation,  where  the  initial 
condition  furnishes  only  one  of  the  required  k+1  values,  or  at  places  where 
the  step  T  is  changed. 

Stability  and  convergence  are  two  important  factors  that  affect  the 
availability  of  a  multiple  step  method.  To  insure  its  stable  and  convergence, 
two  rules  must  be  fulfilled; 

1.   The  characteristic  polynomial 


,k  .  _   Je-1 

+  •  •  •  +  c 
o 


CL  s  +  0L  ,  s"   +•••♦  <x_  -  < 


of  Equation  (3,75)  has  no  root  with  modulus  exceeding  1,  and  the  root  of 
modulus  1  must  be  simple. 
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II.  The  order  of  Che  associated  difference  operator  be  at  least  1. 
A.  Adams -Bashf or th  Mathod 

Consider  the  initial  value  problem 


yf  -  f(t,  y)         y(to)  -  yQ  (3.65) 


An  exact  solution  of  the  differential  equation  (3.65)  by  definition 
satisfies  the  identity 


y(t  +  k)  -  y(t)  - 


f  t-Hc 

f(t,  y)  dt  (3.66) 

t 


for  any  two  values  of  t  in  the  interval  £a,b].  Replace  f (t,y)  in  the  right- 
hand  side  of  Equation  (3,66)  by  an  interpolating  polynomial  on  a  set  of 
points  t  where  y  has  already  been  computed  or  is  just  about  to  be  computed. 
Equate  the  intergral  and  accept  its  value  as  the  increment  of  the  approximate 
values  y  between  t  and  t-tic.   If  it  is  assumed  that  the  interpolating 
points  are  t  ,  t  ,,  t  «>  '*'  tn-a»  ^en   tiie  polynomial  replacing  f (t,y) 
is  given  by 


P(t)  -  L       (-I)*  ("S)vraf  (3.67) 

m-0         ra     p 


c-cp 

s  «  ' 
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where  q    is  an  arbitrary   integer.      This   formula   is  developed    in  Appendix. 
Equation    (3.67)  can  be   substituted    into  Equation    (3.66)    to  yield 


?■■■- 


yp- 


/^p+l 


p(t)dt  - 


r_  vmf  - 


::.-0 


m  v     P 


(3.63) 


m  * 
where  r     -   (-1) 


(    P+l 


("S)ds 
m 


(3.69) 


and  y     -  y(t   ) 

P  P 


Construct  a  generating  function  G(k)  as  follows 


G(k) 


CO 


2    rV 


2     (-1) 

m«0 


t 

f    P+l 


("S)ds 

m 


p+l 


K 


L(-i)ffirs)ds  - 

m 


/Cp+1 


(1  -  k)""  ds 


(3.70) 


The    identity    (l-k)"S  -  exp(-slog(l-k))  causes  Eq.(3.70)   to  become 


G(k) 


(1   -  k)   log(l   -  k) 


this  may  be  written  as 


-  G(k) 


log(l  -  k)       1 


1  -  k 


Since 


1  -  k 


-  1  +  k  +  k  + 
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and 


log(l  -  k) 


1  1  2 

-1  +  -  k  +  -  k  +  ... 

2  3 


one  can  conclude  that 


2  1     *  2 

(r  +  r  k  +  r9k  ■§■•••)(!  +-k+-k  +  •••)»  1+k+k  +  ••• 

2     3 


By  comparing  the  coefficients  of  corresponding  powers  of  k,  a  relation  can 
be  found 


m   0  m-1   _  m-2  ,, 

2       3  m-j-1 


(3.71) 


thus  it  is  possible  to  calculate  r  recursively.   Some  values  of  r 

m  ni 

calculated  from  Eq. (3.71)  are 


m 

0 

1 

2 

3 

4 

5 

6 

r 
m 

1 

1 
2 

5 
12 

3 

8 

251 
720 

95 
288 

19087 

60480 
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If  y^     '    (t)   is  continuous   in  [a,b],  then  y*  (t)  can  be  expressed  as 


y.    .     I     (.!)»(-«)     V<0    +   (-l)(q+1)(-S      )   Tq+1   y(«*2)    (g) 
m-0  m  P  q+1 


where  ^  is  a  point  between  the  largest  and  the  smallest  of  the  values 
t,t  ,  and  t   .   Integrating  yf  between  t  and  t  +1,  we  get 


y(t  ,)  -  y(t  )  -  T   2  r  vmt     +  R 


.AB 


P+l 


;..-0 


m    p    q 


(3.72) 


v/^cre 


AB       (q*D   (q+1) 
I       -  (-1)      T 


f    P+l 


(-s)  y<s+2)  <$>  dt 

q+1         ' 


"s  (q+2) 

since  (   )  is  a  constant  sign  in  the  interval  t  ^t  <t  +1  and  y     (§)  is 

a  continuous  function  of  t,  apply  the  second  mean  value  theorem  of  the 

integral  calculus 


R      -  (-1) 

q 


(q+1)       (q+D       (q+2)        ,         f     P+1      . 


V 


(  '  )  dt 


J  ,.  <1+1 


where  t 


_  _  s  p*  y  t     ..      By  definition  of  r     ,.  this  may  be  written 
p-q  <  ^   <     p+l  J  q+l'  J 
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AB    (q+2)   (q+2) 
R   -  T     y      (?•)  r 
a  7    q+1 


(3.73) 


This  is  the  desired  expression  for  the  remainder  of  the  Adams -Bash forth 
formula. 

B.   Adams -Moulton  tethod. 

This  method  uses  the  form 


P    P-l 


p(t)  dt 


M   *  m 
T  2  r  v  f 

m-0  m 


'p-l 


where 


(3.74) 


*      m 
m 


1   f  P 
T 


y 


Cs)  dt  -  <.»»  1 

m  T 


rs)  ds 

m 


■p-l 


-1 


(3.75) 


the  generating  function  of  the  coefficients  is  determined  as  follows: 


*         *  m 
G  (k)  -  2  r  k 

m-0  m 


log(l  -  k) 


(3.76) 


or 


log(l  -  k)    * 
G  (k)  -  1 


(3.77) 


Expand  the   left-hand  side  of  Eq.    (3.77)   to  power  series 
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*  *  *  it  it 

(l+-k+-k  +•••)(*"     +  r,k  +  r„l<  +  •••    )   «■  1 
2  3  o         1  2 


(3.78) 


It   follows  that 


1     it 


* 


r     +-r     ,+-r     _+•••+  r 

m       ^     m-l        „      m-Z  ,,      o 

2  3  m*l 


1,        in-0 

0,       n-l,2,3,< 


(3.79) 


The  numerical  values  of  r  are  easily  found  from  this  recurrence  relation. 

m 

Some  values  of  r  calculated  from  Eq.  (3.79)  are 


m 

0 

1 

2 

3 

4 

5 

6 

* 

1 

1 

1 

1 

19 

3 

863 

m 

2 

12 

24 

720 

160 

60480 

Since  y  occurs  as  an  argument  in  f  »  f (t  ,  y  )in  the  right  hand 
p  P      P   P 

side  term  of  Eq.(3.74),  it  will  not  be  possible  to  solve  this  equation 
explicitly.  A  better  approach  to  this  solution  is  by  means  of  an  iterative 
procedure. 

Assuming  an  approximation  of  a  solution  of  Eq. (3.74)  has  been  obtained 

to  be  y(o',  calculate  f  ^°'  -  f(t  ,  y^)  and  form  the  difference 
P  P   P 

Vfn01  "  f„   "  f„  i»  "V    f    "  S7f(0)  -  s7f  ,»•'•.   A  better  approxima. 
P      P      P-l       p       p       p-1 


tion  is  then  obtained  from 
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y(1)   -y         +T       I    rvnt™  (3.80) 

P  P-1  m«0     m  p 


(1  \  (1  \ 

Calculating  f«f(t,y       )  and  re -evaluating  the  differences,  a 
P  P       P 


still  better  value  y^2'   is 


(2)  q   *  m  (1) 

y    -y    +T   S  r  V  f  (3.81) 

P      P"1      m»0  m     P 


(r) 

Generally  a  sequence  y   (r  »  0,1,2,3«««)  of  approximations  is 

obtained  recursively  from  the  relation 


y(r)  -y    *T  I    r*vmf(T"l)  (3.82) 

P      P"1      m-0  m  p 


S  ince 


a  (i)     Jn-1   (i)    m-1         (a  "  1»2#,-> 
P  P  P       (i  -  1,2»«0 


(3.83) 


it  follows  that 


m  (r)     id  (r-1)    (r)    (r-1) 
v*  f    -  F  f     -  f    -  f  (3.84) 

P         P       P     P 


Thus 


The  Lipschitz  condition  yields 


Equation    (3.85)   becomes 


or 


where 
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y(r+1>   -  y(r)   -  T     I     r*(^f<r).    /f^l))  (3.85) 

P  P  m-0     m  P  P 


T     2       r      (f  -  f  ) 

n       m       P  P 

m=0 


f(r>  _f(r-l).  .     r        r-1   ,  (3_86) 

P  P  '  '     P        P       ' 


<r*l>         (r)  J      r*  L  ,     <r)  _  y(l>1)  (3>87) 

P  P  m-0       ■  P  P 


Orrt)  .  y<r)    <     (IL  A)r  .     <l)  _     (o)  |  (388) 

P  P  '         P  P 


<5        * 
A  -       £     r 
m-0     m 


The  solution  of  y     of  Eq.  (3.74)    is  now  obtained  by  summing  terms  of 


the   series 


40 


y(0)    +    (ytt)    .   yfe),    +    „(•)    .   /I),    +   ... 


y 
p   p 


*  (y(n)  -  y^)  +...         (3.89) 
P      P 


The  series  on  the  right-hand  side  of  Eq#  (3.89)  will  converge  absolutely 

i    .  ■& 

provided  0  ^  |  TLA   <  1.   In  this  case  the  solution  y  will  exist  and  be 

P 

unique. 

The  local  truncation  error  is 


Rq  -  T«+^  yW'«  (J)  rq+1  (3.90) 


wnere 


p-q  >     p 


C.   Milne's  Method 

Milne's  method  requires  predictor  and  corrector  formulas.   The 
predictor  formula  is  of  the  form 


y     -  y 

nvl    n-3 


t 
(    n+1  4T 

f(t,  y)  dt  =T<2fn_2  -fn-1  +2fn)   (3.91) 

n-3 


and  the  corrector  is  of  the  form 
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y        -  y 
n+:         n-l 


f  fcn+l 


Cn-1 


f (t,  y)  dt  i  -    (f     ,   +  4f     +  f  Al)        (3.92) 

o        n-l  n         n+1 


The   predictor  formula  has  truncation  error 


14       5      (5 ) 
Truncation  error  -  T     y  (£)  (3.93) 

45  ?' 


where 


The  corrector  formula  has  truncation  error 

T5 

Truncation  error  »  y^  '  (f  )  (3.94) 

90        l 


where 

n-l  ^2  **  "n+1 


t    ^§r  ^   t 


The  procedure  is  as  follows: 

Step  1;  Takes  T  small  enough  to  insure  that  the  remainder  term 

ininvolving  is  small  in  the  predictor  formula,  then  find  out  y    from  this 

formula. 

Step  2;  Obtain  a  first  approximation  to  y*   by  substituting  the 

n+1 

value  of  y    obtained  from  step  1  into  the  following  equation 


y?  =  f (t,  y)  y-yatt-t. 

o  o 


42 


Step  3;  Obtain  a  better  approximation  of  y    by  means  of  the  corrector 

n+1 

formula  Eq.  (3.94). 

After  repeating  step  2  and  step  3,  the  value  of  y  .  becomes  very 

accurate.   When  values  of  y  ,  and  f    have  negligible  error,  the  next 

n+1  n+1 

pair  of  values  y     _  and  f     _  may  be  obtained  by  a  repetition  of  the  process. 

let  y^0'   be  the  value  of  y         obtained  from  step  1,  y         be  the  value 
n+1  n+1  n+1 

On) 
of  y         introduced  by  the  corrector  at  first  time,  and  y         be  the  value  of 
n+1  n+1 

y         after  introducing  the  corrector  m  times,  then 
n+1 


(1)  <o)  ,   (1)  (o) 

f         -  f         -  k(y         -  y       )  (3.95) 

n+1         n+1  n+1         n+1 


V7here 


(1)  (o) 

f     .    -  f 

n+x         n+1 


n+1         n+l 


(3.96) 


(o)  (o)  (1)  (1) 

-i_id  f    .,   ra  f<t     .,  y     ,),  f     ,   -  f(t     ,,  y     ,)  (3.97) 

n+1  n+1'  'n+l  n+1  n+1       n+1 


If  f(t,  y)   possesses  a  continuous  first  derivative  f   (t,  y)  with  respect 
to  y,  then  by  the  mean-value  theorem 


lc  -  f(t        ,y       )  (3.98) 

n+1     i  n+1 


wnere 
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<o>  (1) 

;/       lies  between     y       and     y 
(n+1  r.-;-l  n+1 


Equation  (3.92)  together  with  Equation  (3.95)  yield 


(2)  (1)       T        (I)  (o)  (1)  (o) 

n+1         n+1       <*       n+1         n+l  J       L         n+1         n+1   J 


By  the  same  method 


(3)  (2)        ,T        .T       ■     (1)  (o) 

y     ,   -  y     ,    -  (-  k)(-  k)(y         -  y       ) 


where  the  change  of  k   is  negligible. 
Proceeding   in  this  fashion  a  sequence 


n+1  nvl  3  n+1         n+1 


is  obtained. 

Thus 


(o)  (1)  (o)  (2)  (1) 

yA-y+(y-.y-)+(y       -  y     )  +  ••• 

n+.         n+1  n+1         n+1  n+1         n+1 


T  T 


y    :  +  <y       -  y     )  I"  1  +  -  k  +  (-  k)2  +  ...]         (3.100) 

r.v.  n+1  n+1     L  ~  ~  J 


n+1  n+1         n+1 
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T 
Provided  that  T  is  sufficiently  small,  the  value  of  -  k  can.  be  chosen  so 

3 

T 

that  0  ^1"  l<\<£    1  to  insure  convergence  of  the  power  series  on  the  right- 

3 

hand  side  of  Eq. (3.100). 

If  iterations  are  finite,  the  value  of  y  ,  obtained  will  differ  from 

n+1 


y    by 
n*l 


*    ^       (I)    (o)  T   n+1     3 

y   -y   -  (y   -  y   )(-  k>    (  )        (3.101) 

n-M    n+1     n+1    n+1  3        3  -  Tk 


where  n  is  the  number  of  steps  in  the  finite  step  process. 

Let  the  true  values  of  the  y  and  y.8  be  z.  and  z.  respectively.  Define 

iiii 

errors  in  the  y  and  y  by  the  equations 
i      i 


zi  -  y.  +  q%  (3.102) 

z  <=>   y  v  e 
i    i    i 


From  the  differential  equation  (3.64) 


e°   -  z\   -  y]   -  f(t   ,  z   )   -  f(t   ,  y )   -  k   te     -  y  )  (3.103) 

iii  ii  ii  iii 


where     k  -  f (t  ,  *[   )       and      y.<'f,^z 


0 

The  z  and  z  satisfy  the  equation 


45 


5    -  z    -  -  (z*   +  4z'  +  s'   )  -  T^  z(;f;  (  £   }/  90 
r.vi    n-1   3   n-1     n    n+1  -'n+l' 


where 


(3.104) 


t  <  £   <  t 
n   >n+i   n+1 


Subtracting  Equation  (3.94)  from  Equation  (3.104)  results  in 


e    -  e    m   -  (e»    +  4e»  +  e»    )  -  T  A  (3.105) 

/v-.-_    n-1   o   n-1     n    ;v.-l.         n+1 


where 

_<5) 


(^  ,)/90»A  ,  t  <  £  <  t  , 

>n+r        n+1  n  >    n+1 


On  making  use  of  Equation  (3.103)s  this  difference  equation  may  be  written 


1  4  1  5 

(1   -   -  Tk     ,)  e  -   (-  Tk   )  e      +  (1   +  -  Tk      ,)  e      ,    -  T     A     _ 

3       n+1       n+1  0        n'     n  n-1       n-1  n+1 


(3.106) 

Assume    .   k  -  k, 

I  for  i  ■  1,2,«»*  n,  over  a  sufficiently  restricted 

^  A.  -  A 

range. 

Solving  the  difference  equation  yields 


T5A 

e     -CA%  Ci?  + (3.107) 

a  l   l  2   2  <6l  Tk) 

3 
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v;^oro 


2       f~     1     2  2 

-  Tk  +  /l  +  -  T  k  3e 

X  «=  — - -  i  +  Tk  -  e  (3.108) 

i  1  -  -  Tk 

3 


2       /    1  2  2 

-  Tk  -  ,/l  +  -  T  k            Tk       fi 
x     m   J J 3 ~  m   (1 )  S-  .  ey  (3.109) 

2  1  3 

1  -  -  Tk  ' 

3 


1 
G  «  -  Tk  . 

3 


Expressing  the  constants  c„  and  c,.  in  terms  of  e.  and  e  and  setting 

12  1  o 


e     to  be  zero9  Equation   (3.107)  becomes 
o 


n         n         T  A  n 

sA,    -  X_)   +  • (1  -  A.. 


e      £f  -i-  a"  -  C)    +  ' (1   -  C  -   -  9^)  (3.110) 


As  an  approximation  setting  <?1a0!>  Eq,(3.110)  can  be  replaced  by 


5  .  .2 

T  A 


^n-V        ,  ,,n+l  !         "  3  k(tn  '  ^1 
e     »  I  1  -  e      _  +  (-1)  -  Tke     J 


n        2k    ""  2 

(3.111) 


If  k  is  positive,  the  last  term  decreases  as  t  increases.   In  this 

n 

case  the  error  is  always  of  the  same  sign  and  is  multiplied  by  10  each  time 

t  increases  by  2.3/k.   This  case  corresponds  to  a  differential  equation  in 
n 
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which  the  plots  of  the  solutions  for  various  boundary  conditions,  in  the 
neighborhood  of  the  desired  solution,  diverge  to  the  right,  and  the  accumu- 
lated relative  error  may  decrease  even  though  the  accumulated  error  in- 
creases exponentially.   Milne's  method  is  thus  applicable. 

If  k  is  negative,  the  second  term  vanishes,  and  the  error  is  alternately 
positive  and  negative.   This  case  corresponds  to  a  differential  equation 
in  which  the  plots  of  the  solutions  for  various  boundary  conditions  converge 
to  the  right.  The  accumulated  relative  error  increases  rapidly  without 
bound.   This  shows  that  Milne's  method  is  unstable,  and  cannot  be  used. 

Milne's  method  has  two  virtues:   It  supplies  a  running  check  that  the 
method  and  interval  size  are  suitable  and  that  the  computation  is  locally 
accurate  enough  to  warrant  going  on.  Moreover,  only  two  evaluations  of  the 
derivatives  per  step  forward  are  required. 

D.   Hamming  Method 

Milne's  method  is  always  unstable,  because  it  has  two  roots  with  modulus 
one  in  the  characteristic  equation  of  its  corrector  foumula.   Hamming  has 
removed  unstability  by  eliminating  one  of  these  roots  and  thus  modified 
the  corrector  formula  without  changing  the  predictor  formula. 

D.I.  Third  Order  formula 

The  third  order  formula  is  of  the  form 


y    -  ay  +  by    +  T(cy'   +  dy*  +  ey«   )  (3.112) 

n+i     n     n-1      n+1     n     n-1 


Expanding  both  sides  with  respect  to  y  by  Taylor  series  expansion  and 
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2  3 

requiring  exact  fit  for  19T,T  ,T  ,  the  coefficients  a,b9c9d,  and  e  can  be 

determined  to  be 


a  <=>   -4  +  12c  d  =  4  -  8c 

b  «  5  -  12c  e  -  2  -  5c  (3.113) 

c  ™  c 


The  truncation  error  has  the  value 


1  "  3c   4  (4) 
Truncation  error  «  T  y 


The  characteristic  equation  is 


j>2   *   (4  -  12c)/  -  (5  *  12c)  -  0  (3.114) 


Solving  this  equation  yields 


f2a~5+   12c 


Stability  requires 


OSlftU  1 


1  1 
or                 -  >  C>  - 

2  3 
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Setting  c  -  5/12,  Equation  (3.112)  becomes 


7         m   y  +  __  (5y»    +  8y»  -  yt   )  (3.115) 

-v.-_    n   22    n+*     R    "'•'- 


D.  II.  Fourth  Order  Formulas 

This  method  generalises  Milne's  corrector  formula  to  the  form 

y    »  ay  *  by  _  +  cy  _  +  T  [dy«   +  ey«  +  fy'   ]     (3.116) 
n+1     n     n-1     n-2     l      n*l     n     n-1 

and  then  stabilizes  this  formula  by  choosing  suitable  values  of  the 
coefficients  to  minimize  one  of  the  characteristic  roots  with  modulus  one 
in  Milne*s  formula. 

Expanding  Equation  (3.116)  with  respect  y  by  Taylor  series  expansion, 
and  requiring  exact  fit  for  1,  T,  T2,  T3 ,  T4,  yield 


27(1  -  b)  9  -  b 

£   C     ■  '  £   E>   -^— — — • 


24  24 


b  -  b  e  ■  (3.117) 


-3(1  -  b) 

c  -  


e 

18 

+  14b 

24 

-9 

+  17b 

.c 

£ 

24 

24 


In  this  case  the  truncation  error  is  found  to  be 
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5   (5)    "9  *  5b   5   (5) 
kTD  y^  ; r  y^D\  (3.118) 

360 


The  characteristic  equation  of  Equation  (3.116),  using  Equation  (3.117), 
is 

zf   -  9(1  -  b)f2  -  8bf  +  (1  -  b)  -  0.  (3.119) 

The  root  loci  of  f  with  respect  to  b  are  shown  in  Fig.  2.  For  stable 
operation  b  should  lie  in  the  range  -0.6<b<l  to  insure  that  no  charac- 
teristic root  exceeds  one  and  that  the  root  with  modulus  one  is  simple. 
Milne*s  formula  is  a  special  case  of  Equation  (3.116)  with  b-1.   It  is 
easily  seen  from  Fig.  2  that  in  this  case  there  are  two  characteristic  roots 
with  modulus  one. 

For  the  case  b«0,  a  stable  predictor-corrector  formula  can  be  con- 
structed as  follows: 


4T 

predictor       p    -  y  „  +  —  (2y*  -  y«   +  2y»  „) 

n+1    n-3   3    n    n-1     n-2 


112 

modify  m    -  p    -  (p  -  c  ) 

n+1    n+1   12^   n    n 


1 
corrector       c    »  -  (9y  -  y   )  +  3T(f»   -  2y»  -  y*  )     (3.120) 

n+1   8    n    n-2       n+1     n    n-1 


final  value      y    «■  c    +  (p    -  c   ). 

n+1    n+1   121   n+1    n+1 
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Truncation  error  of  corrector  formula  in  this  case  is 


1   5   (5) 
truncation  error  -  -  —  T  y 

40 


as  compared  with  that  of  Milne's  corrector  formula,  there  is  a  125% 
increment. 

In  Hamming  methods  predictor  formulas  are  always  the  same  as  Milne's 
formula  because  if  the  predictor  is  generalized  by  the  same  method  used 
above,  it  would  be  unstable. 

From  the  above  discussion  it  is  easily  seen  that  to  gain  stability  in 
a  predictor-corrector  method  one  must  lose  some  accuracy.   This  loss  in 
accuracy  can  be  compensated  by  shortening  the  interval  of  integration. 
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Fig.    1.      Root  Loci   of  Eq.    (3.114) 

f2  +  (4  -  12c)j°-   (5   +  12c)   -  0 
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Fig.    2.      Root  Loci   of  Eq. (3.119) 


8jJ     -  9(1  -  b)£z  -  8bP  +  (1   -  b)   -  0 
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C.  Comparison  Between  One-Step  Methods  and  Multiple  Step  Methods 

One -step  methods  are  useful  for  obtaining  the  first  few  of  the  solution 
of  a  differential  equation*  but  involve  too  much  labor  to  be  used  for 
obtaining  a  numerical  solution  over  an  extended  range.   However,  multiple 
step  methods  require  special  starting  procedures  which  are  furnished  by 
one-step  methods. 

IV.  Numerical  Transform  Techniques. 

Part  of  the  solution  of  a  linear  differential  equation  with  constant 
coefficients  has  the  form  of  a  convolution.   The  digital  computer  must 
approximate  convolution  on  a  denuraerable  set  of  equally  spaced  points. 
The  tranpezoidal  approximation  is  the  commonly  used  quadrature  method. 

Trapezoidal  approximation  proceeds  as  follows.   Let  g(x)  be  a  function 
defined  and  continuous  in  a  closed  interval  conta inning  x«»0,  and  let  its 
first  four  derivatives  be  continuous  in  the  same  interval,  the  Taylor  series 
representation  of  the  function  with  remainder  is 


2  x3 

g(x)  -g(0)  +xg(1)(0)  +  — g(2)  (0)  +  — g(3)  (0) 

2|  3i 


V 


<*  -^)3    (o) 

-g    (/<)  d/f    (4.1) 


J 


3-' 


Integrating  the  given  function  results  in 
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h 


g(x)  dx 


<g0  +  gx) 


(2) 

o 


h 

12     2 


(2) 


)  +  R(h) 

(4.2) 


The  remainder  R(h)  is  given  by 


h*  r1 


R(h) 


(X  -  2X3  +  X4)  g(4)  (Xh)  dX 


(4.3) 


120 


0<G<1 


A  convolution  can  be  approximated  by  dividing  its  interval  into  a 
finite  number  of  subintervals  each  with  same  width  T  as  follow: 


:  n-1 

fOrt  g(t  -^)  dt»  2 

o  k-o  ; 


(k+l)T 


kT 


f(T)  g(t  -?)   d^   (4.4) 


Trapezoidal  approximation  yields 


f(f)  «(t   -*>  dr-  |    (  T  fk  Sn.k  ♦  fk+1  gn.k.,)        (4.5) 

2        kraU 


Taking  z-transform  of  both  sides  of  Equation  (4.4)  results  in 


■2  [fg  ]  »  T  [  Zf  ][  Zg  ]  -  -  [  f  Zg  +  g  Zf  ] 


(4.6) 
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There  are  many  ways  of  utilizing  trapezoidal  convolution  to  solve  a 
given  differential  equation.   Each  of  these  is  called  a  program.   Three 
classes  of  programs  are  discernible;  The  multiple  integration  substitution 
program,  the  Tustin  integrator  program,  and  the  single  integrator  program. 
Differences  among  them  are  at  the  transitions  from  integration  to  multiple 
integration. 

A.   Tustin  Program 

This  method  solves  an  n-th  order  differential  equation  by  the  state- 
space  spproach,  A  Tustin-like  program  is  demonstrated  as  follows: 
Let  an  n-th  order  differential  equation  of  the  form 


Dny  »  x(t)  (4.7) 


Define  the  row  vectors 


9   /    i       (n-l\ 
w»  -  (y,  y«,  •••  y     ) 


v9  -  (D9  0,  ...  x(t))  (4.8) 


Then 


0      -I 
Dw  +  /  )  w  »  v  (4.9) 


(  )- 

K  0  0  J 


is  equivalent  to  an  nth  order  differentia!  equation.  Take  the  Laplace 
transform  and  divided  by  s  to  obtain 
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111 

w  +  -Aw»-v+-w  (4.10) 

o 
s  s  s 


Employing  Che  sampling  operation  and  trapezoidal    integration  yields 


2 

Zw  +  AZw  -  Aw     »  a2v  -  3v      +  -  8w  (4.11) 

o  o  o 


where 


T       1   +  z  T  /0  -I 

a  ■  -    »         0  " »  a      -     ( 

2       1  -  z  2(1   -  z)  (nxn)  v  0  0 


Since  A  -0f  it  follows  that 


(I  +  aAn)  -i         or  (I  *  aA)  fi  -  (aA)  *  (aA)     +  -«.(-aA) 

-  r  (4.12) 


Then  the  desired  solution  can  be  written  as 


v-    n7.   0   n-1       L     (n-D 

Zy  »  a  Zx  -  £a   x  +  y 

°    1  -  z    o 

2 

T  Z  n   v  o   fn  k^ 

+ 5-  2T  a  "  y        for  n  -  1»2,3"<   (4.13) 

(1  -  zT  k-2      ° 


The  above  equation  is  the  Tustin  program,  if  all  initial  conditions  are 
zero. 

This  program's  great  advantage  is  that  it  finds  the  solution  and  its 
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first  (n-1)  derivatives  in  one  sweep. 

3.  Single -integrator  program. 

Halijak  has  derived  a  single -integrator  program  which  uses  a  sequence 
of  ascending  order  differential  equations  derived  from  the  given  differen- 
tial equation  to  solve  the  same  equation.   Consider  the  differential 

equation. 


n      :\  -, 
D  v  a,  D    *  •  *  •  *  a   ,  D  +  a 
1  n-1     n 


y(t)  ■  x(t)9 


(4. 14) 


First9  S£'c  UP  the  differential  equation 


(D  +  ax)  yL  (t)  -  1         y   (0)  -  0 


(4.15) 


The  Laplace  transform  of  this  differential  equation  yields  after  division 
by  s 


(I  +  _±  )  y  (s) 


(4.16) 


Taking  Z-transform  of  both  sides  and  employing  approximate  trapezoidal 
convolution,  yields 


2Ts 


Zy 


i  Z(- 


1     L1  "z  J[  <2  +  aT>  -  (2  -  -T)  z]        s(s  +  a) 


(4.17) 
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Set  up       r  differential  equation  of  the  form; 


(D   +  axD  +  a2)  y2  (t)  -  1,     y^O)  -  y2(0)  -  0  (4.18) 


taking  Laplace  transform  of  both  sides  and  dividing  by  s(s+a)  yield 


(I    +  )  y(s)  -  — £ 


(4.19) 


s  (s  +  a) 


s  (s  +  a) 


taking  Z-transform  yields 


Zy(s)    +  Ta, 


Z( 


s (s   +  a   ) 


[  Zy(s)]     =    [  Z(-)][Z< 


s(s   +  a) 


•>] 


(A. 20) 


Substituting  equation  (4.17)  into  equation  (4.20)  yields 


T   (1  +  z) 


2Ta 


2y2(s)  -  - 


2   (1  -  z)     (2  +  aT)  -  (4  -  2a  T2)z  +  (2  -  a^z2 


s(s  +  a  s  +  a9) 


(4.21) 


Proceeding   in  this  manner,   let  a  nth  order  differential  equation  be  of  the 

form 


(D  +  a  D  +  •••    *  an_x)  yn>>1    (t)    -  1 


A  recurrence  relation  can  be  constructed  to  be 
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y   (0)  »  y  ,  (0)  =  •••  -  y<n"2)  (o)  -  0  (4.22) 

n-i      n-1  n-1 


T   2  ->  z 
Zy,  (s)  -  -  ( -)  Zy,   (s)  /  [l    +  a  TZy.  ,  (s)J         (4.23) 

1        o    i  1-1       L  1    1-1    J 


2<£i<n 


for  i«n, 


Zy  (s)f  1  +  a  T(Z(y  t  (s))  1  -  T[  Zx(s)  If"  Zy  ,  (s)  ] 

n   L  n     ..->_     J     L       JL   n-1   J 


T 
2 


x  Zy   (s)  } 

o  n-1   J 


(4.24) 


Initial  conditions  other  than  zero  can  be  introduced  at  the  last 
step. 

C.  Multiple-Integrator  Substitution  Program. 

The  multiple  integrator  substitution  program  casts  the  Laplace  trans- 
form of  the  given  differential  equation  into  inverse  powers  of  s  by  a 
suitable  division  and  substitutes  for  them  definite  functions  of  z  defined 
to  oe  e   „ 

C.  I.  Kalijak's  Integrator  Substitution  Program. 

This  method  uses  trapezoidal  convolution  to  generate  z-transform  of 
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(1/s  )  and  (1/s  )f.   The  procedure  is  as  follows: 


for  Z(l/s)  «  1/(1  -  z), 


n  «=  2,        Z(l/sZ)  -  Z(i/s)Z(l/s)T  -  TZ(l/s), 


5/(1  -  2)  , 


(4.25) 


for   n  »  2, 


Z(l/sn)  =  Z(l/sn"1)Z(i/s)T  -  Z(l/sn)T/2, 


T   (1  + 


(— T")  - 


n-1 


1    T 
Z<— s-)  - 

sZ   2 


2   (1  -  z) 


(1  +  z) 


L  a  -  z> 


n-2 


fci  :ion  (4.25)  into  Equation  (4.6)  yields 


1         r      ,  U  T(l  +  2)   1  n"2 

Z(  -  f)  i   T  z/(i  -  z)    I  (Zf  -  0.5f  )   (4.26) 

L   '  U  2(1  -  z)  J 


This  method  yields  moderate  accuracy  approximations  and  is  the  basis 
of  physically  small  computers. 
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SUMMARY 

The  result  of  an  approximate  computation  of  the  solution  of  differen- 
tial equation  is  affected  by  errors.   These  errors  arise  from  different 
causes  and  affect  final  result  in  different  ways. 

Three  types  of  errors  that  are  most  important  are  truncation  error, 
round-off  error,  and  accumulated  error.   Round-off  error  usually  affects 
the  last  retained  digit  of  the  decimal  representation;  its  effect  can  be 
minimized  by  retaining  additional  digits.   Truncation  error  is  due  to 
discarded  terms  in  an  infinite  series.   Sometimes  the  remainder  exceeds  the 
sum  of  terms  retained,  thus  making  the  calculated  result  meaningless.   There- 
-   -. ,  ...  _/.: _..... ce  cf  truncation  error  is  essential.  The  importance  of 
accumulated  errors  depends  on  rate  of  accumulation.   If  the  accumulation 
error  is  unbounded,  the  solution  becomes  meaningless. 

The  search  for  Runge-Kutta  type  formulas  is  important.   It  seems  that 
the  method  given  in  this  report  in  deriving  coefficients  of  Runge-Kutta 
formulas  can  be  applied  to  investigate  six  and  higher  order  formulas. 

Trapezoidal  convolution  using  the  integral  of  the  first  two  terms  of 

Taylor  series  coincides  with  the  average  of  right  and  left  Riemann  sum 

2 
approximations.   The  truncation  error  is  of  order  T  .   Improved  trapezoidal 

convolution  using  higher  order  Taylor  Series  terms  seems  worthy  of  further 

investigation. 

Accuracy  is  not  the  only  consideration  for  evaluating  a  computation 

process.  The  step  size  affects  the  solution  accuracy  as  well  as  the  cost; 

the  more  accurate  solution  is  generally  the  more  expensive.   Frequently  it 

is  desired  to  have  a  solution  within  a  certain  accuracy  with  most  economical 
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computation.   The  selection  of  optimal  size  depends  on  the  method  used  and 
the  problem  solved. 

A  number  of  numerical  transform  techniques  have  been  developed  in  the 
Z-transform  language.   Many  of  these  exhibit  difficulties  for  the  solution 
of  differential  equation  with  non-zero  initial  conditions.   The  programs 
generated  by  Halijak  have  shown  that  proper  re introduction  of  initial  con- 
ditions is  required  for  better  approximation. 
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APPENDIX  1 

Error  Formulas  for  an  Interpolating  Polynomial 

Let  the  function  z (x)  be  defined  on  an  interval  containing  the  q+1 
distinct  points  x,x,«».  x.   It  is  well  known  that  among  all  poly- 
nomials in  x  of  degree  not  exceeding  q  there  exists  exactly  one  polynomial 
P(x)  which  satisfies  the  relations 


P(x  )  -  z(x„)  i  -  0,1,2,.  ..q  (I  -  1) 


The  uniqueness  of  this  interpolating  polynomial  follows  from  the  fact  that 
the  difference  of  any  two  such  polynomials  is  a  polynomial  of  degree   q 
which  has  q+1  zeros  and  therefore  vanishes  identically.  Existence  can  be 
proved  by  exhibiting  the  polynomial  explicitly,  in  the  form 


q      z(x  ) 

?(x)  -  L(x)  2  - (I  -  2) 

i»0  L8(x  )  (x  -  x  ) 


where     L(x)  »  (x  -  x  )<x  -  x, )  ••■•    (x  -  x  ) 

o       1  q 


The  error  committed  in  this  approximation  can  be  estimated  from  the  follow- 
ing Lemmas. 

Lemma  1.   Let  z(x)  have  a  continuous  derivative  of  order  q+1  in  J. 
Then  for  every  point  x  in  J  there  exists  a  point  £  in  the  smallest  interval  1 
containing  both  x  and  the  points  x.  (i«=l,2,«»»q)  such  that 


65 


1  (q+D    > 

z(x)   -   P(x) L(x)  z  H        (J)  (I   -  3) 

(q    +  Di 


Lemma    II.      Let  z(x)   satisfies  the   same  hypothesis  as   in  Lemma   1.      Then 
for  every  x.  (k-0,l,2,««  «q)   there  exists  a  number  £  such  that 


z'(x,  )   -  p»<x,  )   -  z(q+1      (p   L»(x.  )  (I  -  4) 

k  k  (q+D!  '  k 
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There  are  a  number  of  numerical  methods  for  approximating  solution 
of  a  differential  equation.   These  methods  have  two  things  in  common;  the 
calculations  are  performed  with  discrete  values  and  on  a  step-by-step 
basis.  The  purpose  of  this  report  is  to  review  those  methods  that  are  fre- 
quently used  on  digital  computers.  All  these  methods  are  divided  into  two 
classes;  the  classical  numerical  techniques  and  the  numerical  transform 
techniques.   Classical  numerical  techniques  yield  two  types;  one-step  methods 
and  multiple  step  methods.   In  one-step  methods  the  value  of  y  is  solved 
from  its  previous  step.   In  multiple  step  methods  the  value  of  y  is  solved 

from  its  several  previous  steps  y  ,  ,  y  0,  •••  y   •   The  virtue  of  one-step 

n-1*  ^n-2*     •'n-q  K 

methods  is  that  they  are  self  starting  whereas  multiple  step  methods  require 
more  than  one  starting  point  that  they  are  not  self  starting.   However, 
multiple  step  methods  are  more  easier  in  continuing  a  solution  than  one-step 
methods.   The  numerical  transform  techniques  were  developed  recently  by 
Tustin,  Madwed,  Boxer-Thaler,  and  Halijak.   These  methods  use  z-transform 
as  a  language.   Three  classes  are  discernable;  Tustin  program,  single  inte- 
gration program,  and  multiple  integration  substitution  program.   Differences 
among  them  are  the  transition  from  integration  to  multiple  integration. 

Truncation  error  and  accumulated  relative  error  are  important  factors 
that  affect  availability  of  a  method.   It  is  necessary  to  estimate  these 
errors  so  as  not  to  make  computation  meaningless. 


